function ui = cgns_interpolation( grid1,ddc_grid1,ubase1,pbase1, ...
                                  udofs1,pdofs1,uup1, xxx )
% ui = cgns_interpolation( grid1,ddc_grid1,ubase1,pbase1, ...
%                                 udofs1,pdofs1,uup1, xxx )
%
% Interpolate a CG-NS solution from grid1 to the unstructured points
% xxx (as column vectors).
%
% The input arguments grid1, udofs1, pdofs1 and uup1 must be cell
% arrays in order to deal with domain decomposition. Notice that this
% function interpolates also additional fields, such as the vorticity.

 %---------------------------------------------------------------------
 % Set default data
 d = grid1{1}.d;
 xi_default = 1.0/(d+1); % barycenter
 x_default = 0;
 %---------------------------------------------------------------------

 %---------------------------------------------------------------------
 % Compute a smooth vorticity field
 for ig=1:length(grid1)
   omega1{ig} = cgns_vorticity(ubase1,grid1{ig},udofs1{ig},uup1{ig});
 end
 %---------------------------------------------------------------------

 %---------------------------------------------------------------------
 % Interpolate at the points
 igi = 1;
 iei = 1;
 % the assumption here is that the same quad. nodes are used for
 % velocity and pressure
 np = size(xxx,2);
 uue  = zeros(ubase1.pk,d);
 pe   = zeros(pbase1.pk,1);
 ome  = zeros(ubase1.pk,1);
 xi   = zeros(size(xxx,1),1);
 uPHI = zeros(ubase1.pk,1);
 pPHI = zeros(pbase1.pk,1);
 ui.u = zeros(d,np);
 ui.p = zeros(1,np);
 if(d==2)
   ui.w = zeros(1,np);
 end
 tic
 for ip=1:np
   if(mod(ip,floor(np/100))==0)
     toc
     disp(['Point ',num2str(ip),' of ',num2str(np)]);
     tic
   end
   [igi,iei,xii] = locate_point(xxx(:,ip),grid1,ddc_grid1,[igi,iei],1);
   if(iei<0)
     warning(['Unable to locate point ',num2str(xxx(:,ip)'), ...
              '; using default value.']);
     xi(:) = xi_default;
     uue(:,:) = x_default;
     pe(:) = x_default;
     ome(:,:) = x_default;
     igi = 1; iei = 1;
   else
     for id=1:d
       idx = d*(udofs1{igi}.dofs(:,iei)-1) + id;
       uue(:,id) = uup1{igi}(idx);
     end
     idx = d*udofs1{igi}.ndofs + pdofs1{igi}.dofs(:,iei);
     pe(:) = uup1{igi}(idx);
     ome(:) = omega1{igi}( udofs1{igi}.dofs(:,iei) );
     xi(:) = xii(2:d+1);
   end
   % interpolate
   for i=1:ubase1.pk
     uPHI(i) = ev_pol(ubase1.p_s{i},xi);
   end
   for id=1:d
     ui.uu(id,ip) = sum( uue(:,id).*uPHI );
   end
   for i=1:pbase1.pk
     pPHI(i) = ev_pol(pbase1.p_s{i},xi);
   end
   ui.p(ip) = sum( pe.*pPHI );
   % Additional variables:
   % vorticity
   ui.w(ip) = sum( ome.*uPHI );
   % The following would not give a smooth field
   % for i=1:ubase1.pk
   %   for id=1:d
   %     uGPHI(id,i) = ev_pol(ubase1.gradp_s{id,i},xi);
   %   end
   % end
   % coordinate transformation
   % uGPHI = grid1{igi}.e(iei).bi' * uGPHI;
   % if(d==2) % dx(uy)-dy(ux)
   %   ui.w(ip) = uGPHI(1,:)*uue(:,2) - uGPHI(2,:)*uue(:,1);
   % end

 end
 %---------------------------------------------------------------------

return
